#ifndef _CU_FLOAT3_H_
#define _CU_FLOAT3_H_

//**********************************************************************************
// float3 operators and functions
//**********************************************************************************
SDH_CALL float3   operator + (const float3 & u, float v) { float3 r = { u.x + v, u.y + v, u.z + v }; return r; }
SDH_CALL float3   operator - (const float3 & u, float v) { float3 r = { u.x - v, u.y - v, u.z - v }; return r; }
SDH_CALL float3   operator * (const float3 & u, float v) { float3 r = { u.x * v, u.y * v, u.z * v }; return r; }
SDH_CALL float3   operator / (const float3 & u, float v) { float3 r = { u.x / v, u.y / v, u.z / v }; return r; }

SDH_CALL float3   operator + (float v, const float3 & u) { float3 r = { v + u.x, v + u.y, v + u.z }; return r; }
SDH_CALL float3   operator - (float v, const float3 & u) { float3 r = { v - u.x, v - u.y, v - u.z }; return r; }
SDH_CALL float3   operator * (float v, const float3 & u) { float3 r = { v * u.x, v * u.y, v * u.z }; return r; }
SDH_CALL float3   operator / (float v, const float3 & u) { float3 r = { v / u.x, v / u.y, v / u.z }; return r; }

SDH_CALL float3   operator + (const float3 & u, const float3 & v) { float3 r = { u.x + v.x, u.y + v.y, u.z + v.z }; return r; }
SDH_CALL float3   operator - (const float3 & u, const float3 & v) { float3 r = { u.x - v.x, u.y - v.y, u.z - v.z }; return r; }
SDH_CALL float3   operator * (const float3 & u, const float3 & v) { float3 r = { u.x * v.x, u.y * v.y, u.z * v.z }; return r; }
SDH_CALL float3   operator / (const float3 & u, const float3 & v) { float3 r = { u.x / v.x, u.y / v.y, u.z / v.z }; return r; }

SDH_CALL float3   operator - (const float3 & u) { float3 r = { -u.x, -u.y, -u.z }; return r; }

SDH_CALL float3 & operator += (float3 & u, const float3 & v) { u.x += v.x; u.y += v.y; u.z += v.z; return u; }
SDH_CALL float3 & operator -= (float3 & u, const float3 & v) { u.x -= v.x; u.y -= v.y; u.z -= v.z; return u; }
SDH_CALL float3 & operator *= (float3 & u, const float3 & v) { u.x *= v.x; u.y *= v.y; u.z *= v.z; return u; }
SDH_CALL float3 & operator /= (float3 & u, const float3 & v) { u.x /= v.x; u.y /= v.y; u.z /= v.z; return u; }

SDH_CALL float3 & operator += (float3 & u, const float v) { u.x += v; u.y += v; u.z += v; return u; }
SDH_CALL float3 & operator -= (float3 & u, const float v) { u.x -= v; u.y -= v; u.z -= v; return u; }
SDH_CALL float3 & operator *= (float3 & u, const float v) { u.x *= v; u.y *= v; u.z *= v; return u; }
SDH_CALL float3 & operator /= (float3 & u, const float v) { u.x /= v; u.y /= v; u.z /= v; return u; }

SDH_CALL float2 xy         (const float3 & v) { float2 u = { v.x, v.y }; return u; }

SDH_CALL float3 catmullrom (const float3 & P0, const float3 & P1, const float3 & P2, const float3 & P3, float t)
{
  const float ts = t * t;
  const float tc = t * ts;

  return (P0 * (-tc + 2.0f * ts - t) + P1 * (3.0f * tc - 5.0f * ts + 2.0f) + P2 * (-3.0f * tc + 4.0f * ts + t) + P3 * ( tc - ts)) * 0.5f;
}
SDH_CALL float3 lerp      (const float3 & u, const float3 & v, float t ) 
{ 
  return u + t * (v - u); 
}

SDH_CALL float  dot       (const float3 & u, const float3 & v) 
{ 
  return (u.x*v.x + u.y*v.y + u.z*v.z); 
}
SDH_CALL float3 cross     (const float3 & u, const float3 & v) 
{ 
  float3 r = { u.y*v.z-u.z*v.y, u.z*v.x-u.x*v.z, u.x*v.y-u.y*v.x} ; 
  return r; 
}
SDH_CALL float  triple    (const float3 & a, const float3 & b, const float3 & c ) 
{ 
  return dot(a, cross(b, c) ); 
}

SDH_CALL float  length    (const float3 & u) { return sqrt(sqr(u.x) + sqr(u.y) + sqr(u.z)); }
SDH_CALL float3 normalize (const float3 & u) { return u / length(u); }

SDH_CALL float3 clamp     (const float3 & u, float a, float b) 
{ 
  float3 r = { clamp(u.x, a, b), clamp(u.y, a, b), clamp(u.z, a, b) }; return r; 
}
SDH_CALL float3 saturate  (const float3 & u) 
{ 
  return clamp(u, 0.0f, 1.0f); 
}
SDH_CALL float  luminance (const float3 & v) { float3 u = {0.35f, 0.51f, 0.14f}; return dot(u, v); }

SDH_CALL float3 min(const float3 & u, const float3 & v) { float3 r = {min(u.x, v.x), min(u.y, v.y), min(u.z, v.z)}; return r; }
SDH_CALL float3 max(const float3 & u, const float3 & v) { float3 r = {max(u.x, v.x), max(u.y, v.y), max(u.z, v.z)}; return r; }

SDH_CALL float  sad(const float3 & u, const float3 & v) { return (abs(u.x - v.x) + abs(u.y - v.y) + abs(u.z - v.z)); } // Sum of Absolute Differences for xyz components 
SDH_CALL float  ssd(const float3 & u, const float3 & v) { return (sqr(u.x - v.x) + sqr(u.y - v.y) + sqr(u.z - v.z)); } // Sum of Squared  Differences for xyz components 

SDH_CALL float3 set(float3 & v, float x = 0.0f, float y = 0.0f, float z = 0.0f) 
{ 
  v.x = x; v.y = y; v.z = z; return v; 
}
SDH_CALL float3 set(float3 & v, float2 xy, float z = 0.0f) 
{ 
  v.x = xy.x; v.y = xy.y; v.z = z; return v; 
}

SDH_CALL float3 f3()                          { float3 r = {0.0f, 0.0f, 0.0f}; return r; }
SDH_CALL float3 f3(float x)                   { float3 r = {   x,    x,    x}; return r; }
SDH_CALL float3 f3(float x, float y, float z) { float3 r = {   x,    y,    z}; return r; }
SDH_CALL float3 f3(const float2 & u, float z) { float3 r = { u.x,  u.y,    z}; return r; }
SDH_CALL float3 f3(const float4 & u)          { float3 r = { u.x,  u.y,  u.z}; return r; }

SDH_CALL float3 xyZ_to_spherical(const float3 & cartesian)
{
  float r, proj_r, teta, phi; // (r, teta, fi);

  proj_r = sqrtf( sqr(cartesian.x) + sqr(cartesian.y) );
  r      = sqrtf( sqr(cartesian.x) + sqr(cartesian.y) + sqr(cartesian.z) );
  teta   = atanf( cartesian.z / proj_r );
  phi    = acosf( cartesian.x / proj_r ) * sign(cartesian.y);

  float3 spherical = {r, teta, phi};
  return spherical;
}

SDH_CALL float3 spherical_to_xyZ(const float3 & spherical)
{
  float r    = spherical.x;
  float teta = spherical.y;
  float phi  = spherical.z; 

  float x = r * cosf(teta) * cosf(phi);
  float y = r * cosf(teta) * sinf(phi);
  float z = r * sinf(teta);

  float3 cartesian = {x, y, z};
  return cartesian;
}

SDH_CALL float3 xYz_to_spherical(const float3 & cartesian)
{
  float r, proj_r, teta, phi; // (r, teta, fi);

  proj_r =  sqrtf( sqr(cartesian.x) + sqr(cartesian.z) );
  r      =  sqrtf( sqr(cartesian.x) + sqr(cartesian.y) + sqr(cartesian.z) );
  teta   =  atanf( cartesian.y / proj_r );
  phi    =  acosf( cartesian.x / proj_r ) * sign(cartesian.z);

  float3 spherical = {r, teta, phi};
  return spherical;
}

SDH_CALL float3 spherical_to_xYz(const float3 & spherical)
{
  float r    = spherical.x;
  float teta = spherical.y;
  float phi  = spherical.z; 

  float x = r * cosf(teta) * cosf(phi);
  float y = r * sinf(teta);
  float z = r * cosf(teta) * sinf(phi);

  float3 cartesian = {x, y, z};
  return cartesian;
}

#endif